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Abstract. Recent analyses of the helioseismic data have 
produced evidence for a variety of interesting dynamical 
behaviour associated with torsional oscillations. What is 
not so far clear is whether these oscillations extend all 
the way to the bottom of the convection zone and, if so, 
whether the oscillatory behaviour at the top and the bot- 
tom of the convection zone is different. Attempts have 
been made to understand such modes of behaviour within 
the framework of nonlinear dynamo models which include 
the nonlinear action of the Lorentz force of the dynamo 
generated magnetic field on the solar angular velocity. One 
' aspect of these models that remains uncertain is the na- 
£N) ■ ture of the boundary conditions on the magnetic field. 
Here by employing a range of physically plausible bound- 
. ary conditions, we show that for near-critical and mod- 
' erately supercritical dynamo regimes, the oscillations ex- 
£f*^ . tend all the way down to the bottom of the convection 
£S| ' zone. Thus, such penetration is an extremely robust fea- 
ture of the models considered. We also find parameter 
ranges for which the supercritical models show spatiotem- 
poral fragmentation for a range of choices of boundary 
conditions. Given their observational importance, we also 
■^j- ■ make a comparative study of the amplitude of torsional 
CS| ' oscillations as a function of the boundary conditions. 
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1. Introduction 

Recent analyses of the helioseismic data, both from the 
Michelson Doppler Imager (MDI) instrument on board 
the SOHO spacecraft (Howe et al. 2000a) and the Global 
Oscillation Network Group (GONG) project (Antia & 
Basu 2000), have provided strong evidence that the previ- 
ously observed solar torsional oscillations (e.g. Howard & 
LaBonte 1980; Snodgrass, Howard & Webster 1985; Koso- 
vichev & Schou 1997; Schou et al. 1998), with periods of 
about 11 years, penetrate into the convection zone (CZ) 
to depths of at least 10 percent in radius. 

These studies have also produced rather conflicting re- 
sults concerning the dynamical behaviour near the bottom 
of the convection zone. Thus Howe et al. (2000b) find ev- 
idence for the presence of torsional oscillations near the 
tachocline situated close to the bottom of the convection 
zone, but with a markedly shorter period of about 1.3 
years, whereas Antia & Basu (2000) do not find such os- 
cillations. Given the uncertainties in the helioseismic data, 
what is not certain so far is (i) whether torsional oscilla- 
tions do extend all the way to the bottom of the CZ and 
(ii) whether there are different oscillatory modes of be- 
haviour at the top and the bottom of the CZ. 

Work is in progress by a number of groups to re- 
peat these analyses in order to answer these observational 
questions. In parallel, attempts have been made to ap- 
proach these questions theoretically by modelling varia- 
tions in the CZ within the framework of nonlinear dynamo 
models which include a nonlinear action of the azimuthal 
component of the Lorentz force of the dynamo generated 
magnetic field on the solar angular velocity (Covas et al. 
2000a, b; Covas et al. 2001a, b, see also erratum in Covas et 
al. 2002). According to these results, for most ranges of dy- 
namo parameters, such as the dynamo and Prantdl num- 
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bers, the torsional oscillations extend all the way down to 
the bottom of the convection zone. In addition, spatiotem- 
poral fragmentation/bifurcation (STF) has been proposed 
as a dynamical mechanism to account for the possible ex- 
istence of multi-mode behaviour in different parts of the 
solar CZ (Covas et al. 2000b, 2001a,b, 2002). In all these 
studies the underlying zero order angular velocity was cho- 
sen to be consistent with the recent helioseismic data. 

As in much astrophysical modelling, an important 
source of uncertainty in these models is the nature of 
their boundary conditions. Given this uncertainty, and the 
fact that boundary conditions can alter qualitatively the 
behaviour of dynamical systems, it is important to see 
whether employing different boundary conditions can sig- 
nificantly change the dynamics in the CZ, and in particu- 
lar whether the two dynamical modes of behaviour men- 
tioned above are robust with respect to plausible changes 
in the boundary conditions. This is important for two rea- 
sons. Firstly, in the absence of precise knowledge about 
such boundary conditions, it is important that the dy- 
namical phenomena of interest predicted by such mod- 
els can survive reasonable changes in ill-known boundary 
conditions. Secondly, it may in principle be possible for 
qualitative changes found as the boundary conditions are 
altered to be used as a diagnostic tool to determine the 
range of physically reasonable boundary conditions in the 
solar context. 

Here, by considering a number of families of boundary 
conditions, we show that the penetration of torsional os- 
cillations to the bottom of the CZ is indeed robust with 
respect to a number of plausible changes to the bound- 
ary conditions. We also find spatiotemporal fragmenta- 
tion in these models with a variety of, but not all, choices 
of boundary conditions. Given the observational impor- 
tance of the amplitudes of the torsional oscillations, we 
also make a comparative study of their magnitudes as a 
function of the boundary conditions. 
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2. The model 

We shall assume that the gross features of the large scale 
solar magnetic field can be described by a mean field dy- 
namo model, with the standard equation 
dB 

= V x (u x B + aB - 77V x B). (1) 



iV?7, the term proportional to V77 repre- 



dt 

Here u 

sents the effects of turbulent diamagnetism, and the ve- 
locity field is taken to be of the form v = v + v', where 
vq = QorsmO, Slo is a prescribed underlying rotation law 
and the component v' satisfies 
dv' _ (V x B) x B 
~dt~ . 



.0 + vD 2 



(2) 



Ho pr sm 

where D 2 is the operator ^ + \ -r ^- sin e , 
■^-g) and /io is the induction constant. The assumption 
of axisymmetry allows the field B to be split simply into 
toroidal and poloidal parts, B = Br+Bp = B(f>+\7 xA(f>, 
and Eq. ([!]) then yields two scalar equations for A and 
B. Nondimensionalizing in terms of the solar radius R 
and time R 2 /i]o, where 770 is the maximum value of 77, 
and putting fl = a = cvotv, 77 = 770^?, B = BqB 

and v' = Q*Rv', results in a system of equations for 
A, B and v'. The dynamo parameters are R a = ao-R/770, 
R u = fl*R 2 /770, P r — vo/tiq, and 77 = 77/770, where il* is 
the solar surface equatorial angular velocity. Here vq and 
770 are the turbulent magnetic diffusivity and viscosity re- 
spectively and P r is the turbulent Prandtl number. Our 
computational procedure is to adjust R u so as to make 
the cycle period be near the solar cycle period of about 
22 years for the marginal dynamo number, and then to 
allow R a and, to some extent, P r to vary. The density p 
is assumed to be uniform. 

Eqs. ([[]) and (g) were solved using the code described 
in Moss & Brooke (2000) (see also Covas et al. 2000b) to- 
gether with the boundary conditions given below, over the 
range ro < r < 1, < 6 < ir. We set ro = 0.64, and with 
the solar CZ proper being thought to occupy the region 
r > 0.7, the region ro < r < 0.7 can be thought of as an 
overshoot region/tachocline. In the following simulations 
we used a mesh resolution of 61 x 101 points, uniformly 
distributed in radius and latitude respectively. 

In this investigation, we took f^o in 0.64 < r < 1 to be 
given by an interpolation on the MDI data obtained from 
1996 to 1999 (Howe et al. 2000a). We set a = a r (r)f(6), 
where f{0) was chosen to be sin 2 cos 9 or sin 4 9 cos 9. The 
angular structure of a is quite uncertain, and both these 
forms have been used in the literature (see e.g. Riidiger & 
Brandenburg 1995) and their choice here is simply to make 
the butterfly diagrams more realistic. We took a r = 1 in 
all or part of the CZ (see below for details), with cubic 
interpolation to zero at r = 7*0 and r = 1 in the cases 
where a r ^ 1 everywhere. Throughout we take a r > 
and R a < 0. Also, in order to take some account of the 
likely decrease in the turbulent diffusion coefficient 77 in 



the overshoot region, we allowed a simple linear decrease 
from 77 = 1 at r — 0.8 to 77 = 0.5 in r < 0.7. 



3. The choice of boundary conditions 

Boundary conditions on magnetic fields are often rather 
ill-determined when modelling astrophysical systems. 
This is certainly true in the case of the Sun and solar-type 
stars. Given this uncertainty, we shall consider a number 
of physically motivated families of boundary conditions 
and investigate the consequences of each on the dynamics 
of the CZ. In particular we shall study whether they al- 
low penetration of torsional oscillations all the way to the 
bottom of the CZ as well as supporting spatiotemporal 
fragmentation. We note that at 9 = and ir symmetry 
conditions imply A = B = 0. In this article we shall con- 
centrate on the changes to the outer boundary conditions 
only. 



3.1. Boundary conditions at r = 7*0 

The detailed physics is uncertain near the base of the 
computational region (r = ro). Given that the angular 
momentum flux out of a region with boundary S from 
the magnetic stresses is J s (BBrsin9)dS, we set B = 
on r = ro in order to ensure zero angular momentum 
flux across the boundary and, correspondingly, stress-free 
conditions were used for v'. The condition dA/dr = A/6 
crudely models A falling to zero at distance 5 below r = ro 
(cf. Moss, Mestel & Tayler 1990; Tworkowski et al. 1998). 
We chose 6 = 0.03, but the general nature of the results 
is insensitive to this choice. Taking 6 > is computation- 
ally helpful as it reduces somewhat the field gradients near 
r = ro, although it is not essential. 



3.2. The outer boundary conditions 

At the outer boundary r = R, we shall, in view of the 
uncertainties regarding the outer boundary conditions, 
consider a number of different but physically reasonable 
choices. 

One of the common choices for the outer boundary 
conditions adopted in literature is the 'vacuum' bound- 
ary condition, in which the poloidal field within r = R is 
smoothly joined, by a matrix multiplication, to an exter- 
nal vacuum solution; the azimuthal field B = 0. 

Given the dynamic nature of the solar surface, the 
usual vacuum conditions can, to some extent at least, 
be regarded as a mathematically convenient idealization. 
Some aspects of this issue have recently been discussed at 
length by Kitchatinov, Mazur & Jardine (2000), who de- 
rive 'non-vacuum' boundary conditions on both B and 
Bp. We also consider families of boundary conditions 
which deviate from the vacuum conditions and refer to 
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these as 'open'. As a convenient and flexible general form 
for the boundary conditions at the surface, we write 

dA „ dB 

r— + ni A = 0, r—+n 2 B = 3 

dr dr 

where n\ and n 2 are constants that parameterize the 

boundary conditions and, to some extent, their degree of 
openness. With n\ — 1, n 2 = the two conditions for 
A and B reduce to ^£r^ = Bg = and 4^ = respec- 
tively. The condition Bg = has been adopted previously 
by some investigators. The limit n 2 — > oo gives the of- 
ten used B — 0. As n± increases, the penetration of the 
poloidal field through the surface decreases, and in the 
limit m — > oo the boundary condition is A = 0, and all the 
poloidal field lines then close beneath the surface r = R, 
which is thus the limiting field line. Using the vacuum 
boundary condition for Bp gives poloidal field lines that 
mostly make a modest angle with the radial direction, and 
so we can anticipate that, by taking small values of m and 
large values of n 2 , we will obtain solutions that resemble 
in some ways those found by using the vacuum boundary 
conditions mentioned above. (But note that Eq. (|^) gives 
strictly local conditions on the field components, whereas 
the vacuum condition on the poloidal field is essentially 
nonlocal.) 

We note one further technical point. The angular mo- 
mentum flux through r = R is non-zero if both Bp and B 
are non-zero there. Whilst the Sun certainly is losing an- 
gular momentum, we are not trying to model this process 
here, and so will only consider models in which the angular 
momentum 'drift' of the dynamo region is small enough 
that we can consider it to be a unchanging background for 
the dynamo calculations. 

Now in order to find the range of values of n\ and n 2 
such that the resulting 'partially open' boundary condi- 
tions are physically plausible, we need to ensure that the 
chosen boundary conditions result in appreciable poloidal 
flux penetrating the surface. Thus we calculated the aver- 
age over the dynamo cycle of the ratio of the flux of the 
poloidal field at the surface to the corresponding value 
within the CZ, given by 

_ [Rsm(6)A(R,G)] 

Surface / A \ 

s= [r S in(9)A(r,9)] lnncr ' U 
as a function of n\ say. Here R is the model radius, the 
numerator is evaluated at the surface r — R, and the 
denominator is evaluated inside the dynamo region ('CZ') 
tq < r < R. As n\ increases we would expect this ratio 
to decrease. We consider a boundary condition to be, in 
principle, viable if the ratio F s is not too small compared 
with the corresponding value in the 'standard' vacuum 
case. 

4. Results 

Using the above model, we studied the dynamics in the 
convection zone subject to three sets of boundary condi- 
tions; namely, the vacuum boundary condition and two 
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Fig. 1. The radial (r-t) contours of the angular veloc- 
ity residuals 5£l as a function of time for a cut at 25 
degrees latitude, with vacuum boundary conditions and 
R a = _5.5, p r = 0.6, R u = 60000. Note the fragmenta- 
tion at the bottom of the convection zone and the resulting 
difference in periods of oscillations at the top and at the 
bottom. Darker and lighter regions represent positive and 
negative deviations from the time averaged background 
rotation rate. 



families of open boundary conditions which we shall refer 
to as boundary conditions (1) and (2). Also in order to 
demonstrate that penetration of the torsional oscillations, 
as well as spatial fragmentation, can occur with various 
changes in other ingredients of the model, we have chosen 
examples with different forms for these ingredients. 

4-1- Vacuum boundary conditions 

With this choice of boundary conditions, we found that 
for critical and moderately supercritical regimes, the tor- 
sional oscillations extend all the way down to the bottom 
of the CZ (see also Covas et al. 2000a, where the critical 
value of the dynamo number for the onset of dynamo was 
found to be R a k, —3.16). In addition we found ranges of 
dynamo parameters for which supercritical models showed 
spatiotemporal fragmentation. As an example of such frag- 
mentation, we show in Fig. [l]the radial contours of the an- 
gular velocity residuals Sil as a function of time for a cut at 
25 degrees latitude. In this case we took f(9) = sin 4 9 cos 9, 
with a r = 1 throughout the computational region. The 
parameter values used were R a = —5.5, P r = 0.6 and 
Ruj — 60000. We also show in Fig. || the magnetic butter- 
fly diagram. 

We note that this is the first time STF has been ob- 
tained with vacuum boundary conditions. This is despite 
the fact that the range of parameter values for which STF 
is present for this model is rather wide, as is the case with 
the boundary conditions considered below. What seems 
to occur here is that the onset of spatiotemporal fragmen- 
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Fig. 2. Butterfly diagram of the toroidal component of 
the magnetic field B at tq — 0.90i?, with the vacuum 
boundary conditions and the parameter values given by 
R a = _5.5, P r = 0.6, R u = 60000. Dark and light shades 
correspond to positive and negative values respectively. 



tation is close to a bifurcation point, which disrupts the 
butterfly diagram and rather confuses the diagnosis of the 
situation. 

4-2. Open boundary conditions (1) 

We now take boundary conditions at r = R to be given 
by Eq. (|j|), where n-i = and n\ is varied. 

In order to find the range of values of n\ which can be 
considered as physically plausible, we calculated the ratio 
of the poloidal fluxes F s given by (Q) as a function of n±. 
The results are given in Fig. ||. To give an idea of how this 
ratio compares to that of models with vacuum boundary 
conditions, we also calculated F s with vacuum boundary 
conditions with the same parameter values as for Fig. |l|, 
and found that F s ~ 0.25. As can be seen from Fig. [3|this 
is comparable in magnitude to the values obtained with 
the above open boundary conditions for a wide range of 
values of n\ given by n\ < 40. 

With these open boundary conditions, we again found 
that for slightly and moderately supercritical dynamo 
regimes, the torsional oscillations extend all the way down 
to the bottom of the CZ. There are also ranges of dy- 
namo parameters for which spatiotcmporal fragmentation 
is found. When n\ is close to 1 there is appreciable an- 
gular momentum drift, and so we do not consider these 
solutions further here. As an example of a case with STF 
(and with negligible angular momentum drift over the in- 
tegration interval), we show in Fig. ^ the radial contours of 
the angular velocity residuals 6£l as a function of time for 
a cut at 25 degrees latitude. In this case f(8) = sin 4 9 cos 9 
and a r is given by a r = 1 for 0.7 < r < 0.8 with cubic 
interpolation to zero at r = tq and r = 1. The parameter 
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Fig. 3. F s as a function of n\ for the open boundary 
conditions (1) given by equation (H), with n2 = 0. 



1 .00 



0'j 

D 
01 



0.90 



0.80 



0..70 




10 20 30 
Time (years) 



40 



Fig. 4. The radial (r-t) contours of the angular velocity 
residuals <5f2 as a function of time for a cut at 25 degrees 
latitude, with the boundary conditions given by (|) with 
ni = 25, R a = -5.0, P r = 0.7, i? w = 50000. Note the 
fragmentation at the bottom of the convection zone and 
the resulting difference in periods of oscillations at the top 
and at the bottom. Darker and lighter regions represent 
positive and negative deviations from the time averaged 
background rotation rate. 



values used are R a = —5.0, P r = 0.7, R u = 50000 and 
m = 25. 

Fig. U shows the angular velocity residuals at R = 0.92 
(i.e. the near-surface torsional oscillations), with latitude 
and time, for the same parameter values as in Fig. |^. The 
poloidal field lines and toroidal field contours for this case 
are presented in Figs. ^ and ^. In Fig. |7], the effect of 
the open boundary condition on B is seen to be relatively 
minor - most of the toroidal field is concentrated near the 
tachocline, and far from the surface. The poloidal field 
(Fig. ^) is more uniformly distributed. 
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Fig. 5. Angular velocity residuals at R = 0.92 with lati- 
tude and time, with boundary conditions given by Eq. (||), 
m = 25, R a = -5.0, P r = 0.7, i? w = 50000. A tempo- 
ral average has been subtracted to reveal the migrating 
banded zonal flows. Darker and lighter regions represent 
positive and negative deviations from the time averaged 
background rotation rate. 




Fig. 6. The poloidal field lines; same parameters as in 
Fig. §. 



4-3. Open boundary conditions (2) 

We now consider boundary conditions given by Eq. (||), 
with n\ of order one and 1 < n 2 < 500 (so B « for large 
n 2 ). We found that for m values close to one our model 
again has significant angular momentum drift, which in- 
creases as n\ decreases to 1. For larger values of ri\ the 
drift is negligible, and the ratio F s is again 'reasonable', 
being comparable with the vacuum case. We found that 
putting rii = 2 gave satisfactory behaviour, and this is the 
case that we discuss in detail below. 



Fig. 7. The toroidal field contours Bp; same parameters 
as for Fig. |[ 
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Fig. 8. The radial (r-t) contours of the angular velocity 
residuals 5£l as a function of time for a cut at 20 degrees 
latitude, with the open boundary conditions (2) and R a = 
-16, P r = 0.22, R u = 48000, m = 2 and n 2 = 400. Note 
the fragmentation at the bottom of the convection zone 
and the resulting difference in periods of oscillations at 
the top and a region near the bottom of the CZ around 
ro = 0.7. Darker and lighter regions represent positive and 
negative deviations from the time averaged background 
rotation rate. 



For all such cases we again found that for slightly and 
moderately supercritical dynamo numbers, the oscillations 
extend all the way down to the bottom of the CZ. In ad- 
dition we found ranges of dynamo parameters for which 
the supercritical model show spatiotemporal fragmenta- 
tion. An example, which has negligible angular momen- 
tum drift over the integration interval, is shown in Fig. ^| 
Here a r = 1 for 0.7 < r < 0.8 with cubic interpolation 
to zero at r = tq and r — 1 and f{9) = sin 2 flcosf?. The 
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parameter values used were R a = — 16.0, P r = 0.22 and 
R u = 48000, and boundary conditions were given by (||) 
with n\ = 2 and n 2 = 400. 

4-4- Amplitudes of oscillations as a function of boundary 
conditions 

Another important issue from an observational point of 
view is the way the amplitudes of the torsional oscillations 
vary as a function of model ingredients and parameters, 
as well as with depth in the CZ. 

To begin with, we verified that for a given boundary 
condition the amplitudes of the oscillations increase as 
the dynamo number R a and the Prantdl number P r are 
increased (see also Covas et al 2001a). For orientation, we 
recall that the observed surface amplitudes in the case of 
the Sun are latitude dependent and of order of one nHz 
(see e.g. Howe et al. (2000b)). 

We also made a comparative study of the amplitudes as 
a function of changes in the boundary conditions. Briefly 
we found that typically the solutions with vacuum bound- 
ary conditions have rather smaller amplitudes of oscilla- 
tions, especially near the surface. For example, for the 
model of Fig. 1 (which has a high R a , spatiotcmporal 
fragmentation and thus would be expected to have higher 
amplitudes), we found the mean averaged amplitudes to 
be 0.72,0.19,0.09 nHz at the depths r = 0.70,0.88 and 
0.95 respectively. 

For the models with open boundary conditions, we 
found the amplitudes to be on average higher than the 
vacuum case, specially near the surface. We have sum- 
marised in Fig. H our calculations of the amplitudes of 
oscillations for the models with open boundary conditions 
given by Eq. (||) and n 2 = 0, for a range of n\ given 
by 1 < n\ < 150. Here the dynamo parameters were 
R a = -5.0, P r = 0.7, R u = 50000 and f {9) = sin 4 6>cos6>, 
with a r — 1 for 0.7 < r < 0.8, with cubic interpolation to 
zero at r = tq and r = 1. As can be seen, the amplitudes 
grow at all depths with increasing n\ and saturate around 
n\ ~ 50. As an example, the models with n\ values around 
ni = 25, (corresponding to Fig. ||, with STF), have ampli- 
tudes that are more than double those found above with 
vacuum boundary conditions, down to the level ro = 0.88. 

5. Discussion 

We have made a detailed study of the effects of the bound- 
ary conditions on the dynamics in the solar convection 
zone, by employing various forms of outer boundary con- 
ditions. 

In all the models considered here (as well as other re- 
sults not reported), we find that in near-critical and mod- 
erately supercritical dynamo regimes the torsional oscilla- 
tions extend all the way down to the bottom of the CZ. In 
this way our results, taken altogether, demonstrate that 
such penetration is extremely robust with respect to all 
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Fig. 9. The mean amplitudes of torsional oscillations, at 
radii ro — 0.70, 0.88 and 0.95. in models with open bound- 
ary conditions given by Eq. (|^) and n 2 — 0, with n\ rang- 
ing from 3 to 150. 



the changes we have considered both to the boundary 
conditions, and the dynamo parameters such as the dy- 
namo and Prantdl numbers, in addition to variations in 
the model ingredients such as the a and n profiles and the 
rotation inversion. 

We deduce, that if our dynamo model (which is ba- 
sically a standard mean field dynamo) has any validity, 
then observers should expect to find that the solar tor- 
sional oscillations penetrate to the tachocline. However, 
given the significant uncertainties that still remain in he- 
lioseismic measurements, especially the limited temporal 
extent of the data available, this issue may not be resolv- 
able at present (see e.g. Vorontsov et al. 2002). 

In all cases we have found supercritical dynamo 
regimes with spatiotemporal fragmentation for a range of, 
but not all, dynamo parameters. This, together with our 
previous work, shows that fragmentation occurs with a 
variety of forms of a (and also that it is not confined to 
a particular inversion for the solar angular velocity). For 
still more supercritical dynamo regimes we find a series of 
spatiotemporal fragmentations, leading eventually to spa- 
tiotemporal chaos, i.e. disappearance of coherence in the 
dynamo regime. 

These results are of potential importance in interpret- 
ing the current observations, especially given their diffi- 
culty in resolving the dynamical regimes near the bottom 
of the convection zone. However given the variety of dy- 
namical behaviour possible theoretically near the bottom 
of the CZ, we cannot comment definitively on the reported 
1 — 3 yr oscillations. 

Finally, given the observational importance of the am- 
plitudes of the torsional oscillations, we have made a com- 
parative study of their magnitudes as a function of the 
boundary conditions. We found that on average the am- 
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plitudes are smaller for the models with vacuum boundary 
conditions than for those with open boundary conditions. 
An important ingredient that our model omits and which 
seems bound to have an effect on the amplitudes of tor- 
sional oscillations throughout the CZ, is that of density 
stratification. We intend to return to this issue in a future 
publication. 
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